%%
%% Copyright 2007, 2008, 2009 Elsevier Ltd
%%
%% This file is part of the 'Elsarticle Bundle'.
%% ---------------------------------------------
%%
%% It may be distributed under the conditions of the LaTeX Project Public
%% License, either version 1.2 of this license or (at your option) any
%% later version.  The latest version of this license is in
%%    http://www.latex-project.org/lppl.txt
%% and version 1.2 or later is part of all distributions of LaTeX
%% version 1999/12/01 or later.
%%
%% The list of all files belonging to the 'Elsarticle Bundle' is
%% given in the file `manifest.txt'.
%%

%% Template article for Elsevier's document class `elsarticle'
%% with harvard style bibliographic references
%% SP 2008/03/01
%%
%%
%%
%% $Id: elsarticle-template-harv.tex 4 2009-10-24 08:22:58Z rishi $
%%
%%
\documentclass[preprint,authoryear,12pt]{elsarticle}

%% Use the option review to obtain double line spacing
%% \documentclass[authoryear,preprint,review,12pt]{elsarticle}

%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%% \documentclass[final,authoryear,1p,times]{elsarticle}
%% \documentclass[final,authoryear,1p,times,twocolumn]{elsarticle}
%% \documentclass[final,authoryear,3p,times]{elsarticle}
%% \documentclass[final,authoryear,3p,times,twocolumn]{elsarticle}
%% \documentclass[final,authoryear,5p,times]{elsarticle}
%% \documentclass[final,authoryear,5p,times,twocolumn]{elsarticle}

%% if you use PostScript figures in your article
%% use the graphics package for simple commands
%% \usepackage{graphics}
%% or use the graphicx package for more complicated commands
%% \usepackage{graphicx}
%% or use the epsfig package if you prefer to use the old commands
%% \usepackage{epsfig}

%% The amssymb package provides various useful mathematical symbols
%\usepackage{amssymb}
%% The amsthm package provides extended theorem environments
\usepackage{amsthm}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{colortbl}

%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers after \end{frontmatter}.
%% \usepackage{lineno}

%% natbib.sty is loaded by default. However, natbib options can be
%% provided with \biboptions{...} command. Following options are
%% valid:

%%   round  -  round parentheses are used (default)
%%   square -  square brackets are used   [option]
%%   curly  -  curly braces are used      {option}
%%   angle  -  angle brackets are used    <option>
%%   semicolon  -  multiple citations separated by semi-colon (default)
%%   colon  - same as semicolon, an earlier confusion
%%   comma  -  separated by comma
%%   authoryear - selects author-year citations (default)
%%   numbers-  selects numerical citations
%%   super  -  numerical citations as superscripts
%%   sort   -  sorts multiple citations according to order in ref. list
%%   sort&compress   -  like sort, but also compresses numerical citations
%%   compress - compresses without sorting
%%   longnamesfirst  -  makes first citation full author list
%%
%% \biboptions{longnamesfirst,comma}

% \biboptions{}

\journal{Computers and Mathematics with Applications}

\begin{document}

\begin{frontmatter}

%% Title, authors and addresses

%% use the tnoteref command within \title for footnotes;
%% use the tnotetext command for the associated footnote;
%% use the fnref command within \author or \address for footnotes;
%% use the fntext command for the associated footnote;
%% use the corref command within \author for corresponding author footnotes;
%% use the cortext command for the associated footnote;
%% use the ead command for the email address,
%% and the form \ead[url] for the home page:
%%
%% \title{Title\tnoteref{label1}}
%% \tnotetext[label1]{}
%% \author{Name\corref{cor1}\fnref{label2}}
%% \ead{email address}
%% \ead[url]{home page}
%% \fntext[label2]{}
%% \cortext[cor1]{}
%% \address{Address\fnref{label3}}
%% \fntext[label3]{}

\title{A comparison study of Multi-Component Lattice Boltzmann models}

%% use optional labels to link authors explicitly to addresses:
%% \author[label1,label2]{<author name>}
%% \address[label1]{<address>}
%% \address[label2]{<address>}

\author[jy]{Jianhui Yang}
\ead{jianhui.yang10@imperial.ac.uk}
\author[jy]{Edo Boek\corref{cor1}}%\fnref{fn1}
\ead{e.boek@imperial.ac.uk}

\cortext[cor1]{Corresponding author}
%\fntext[fn1]{This is the specimen author footnote.}
\address[jy]{Department of Chemical Engineering, Imperial College London, London, SW7 2AZ, United Kingdom}


\begin{abstract}
A comparison study of three different multi-component Lattice Boltzmann models is carried out to explore their capability of describing binary immiscible fluid systems. The Shan-Chen pseudo potential model, the Oxford free energy model, and the Colour gradient model are investigated using the multi-relaxation time scheme (MRT) algorithm to study the flow of oil, water and supercritical CO2 in reservoir rocks. Poiseuille flow of layered immiscible binary fluids and capillary fingering phenomena were examed and evaluated against analytical solutions. In addition, we examine the capability of the various models to simulate fluids with significant viscosity and density contrast and suitable interface thickness.The Shan-Chen model is found capable of simulating high density ratio fluids, but gives relatively low numerical stability and thick interfaces. Both the Colour gradient model and Oxford free energy model are capable of simulating viscosity contrast fluids, and can recover the analytical solutions of Poiseuille flow, capillary rise and fingering simulations. The Oxford free energy model can predict any diffusion-convection binary fluids system with solid thermodynamical foundations. The physical foundations of the Colour gradient model are less clear so that the generalisation to other  compositions is less straightforward.We conclude that the Oxford free energy model seems to be the most appropriate model to simulate the flow of water, supercritical CO2 and oil in porous media, due to its capability of simulating binary fluids with high viscosity contrast and  high numerical stability.


\end{abstract}

\begin{keyword}
%% keywords here, in the form: keyword \sep keyword
Lattice Boltzmann method, multi cpmponent models, binary fluids
%% MSC codes here, in the form: \MSC code \sep code
%% or \MSC[2008] code \sep code (2000 is the default)

\end{keyword}

\end{frontmatter}

% \linenumbers

%% main text
\section{}
\label{}

\section{Introduction}
\subsection{Background}
Multi-component flows are of great importance in many engineering fields including petroleum industrial, biochemical engeering, chemical engineering. Several conventional CFD techniques including volume of fluid (VOF) method, and the level-set method have been used to study the multi-component flow. Interface dynamics at large scale can be captured by these techniques, but a huge number of information of tiny interfaces is missing \cite{guo_book}.The lattice Boltzmann method can be an alternative solution for simulations of complex flow due to its statistical physics background, ease with implementation, strong strength of dealing with complex geometries and  inherent parallelism.\\

\subsection{Multi-component lattice Boltzmann models}
The lattice Boltzmann multi-compnent medel \cite{doolen_1990} \cite{chen_1992} was firstly introduced to multi-component flow simulation by Gunstensen et al (1991) \cite{color1991} with the name colour gradient model.Two components represent two types of fluid with their own distribution functions, and follow their own evolution equation. They are named as red particles and blue particles. The collision step includes self-interactions and cross-interactions with other types of particles. A colour function gradient was introduced to calculate the surface tension between different phases \cite{grunau_1993}. To segregate the phases, mixing near the interface should be minimized. A procedure called recoloring is proposed for this minimizing process \cite{tolke_2002} \cite{ortona_1995}. \\


The free energy model was developed by Swift et al (1996) \cite{freeenergy}. This model includes thermodynamic equilibrium functions of phases, and a term describing the surface tension is added to the equilibrium distribution function. This allows the free energy model to specify the surface tension more easily than other multiphase multi-component models.And it is a fully thermodynamically consistent binary fluid lattice Boltzmann model.To reduce the spurious currents, Pooley et al (2008) \cite{pooley_2008} proposed a modified distribution functions for the Free Energy Model that the spurious velocities may be decreased by an order of magnitude compared to previous models.\\

A pseudo potential lattice Boltzmann model was developed by Shan and Chen (1992) \cite{shanchen} \cite{shanchen1995}. The principal characteristic of this model is an interaction force between particles is introduced to have a consistent treatment of the equation of state for a non-ideal gas.Shan (2006)reported that the spurious velocity is due to the lack of sufficient isotropy in the calculation of the gradient term for the interaction force. Finite difference gradient operators with higher order of isotropy was proposed and the spurious current was found decreasing significantly. \cite{shan2006}\\


\section{Comparison of Multi-component lattice Boltzmann models}
We briefly review multi-component lattice Boltzmann models including the Shan-Chen pseudo potential model, the Free Energy model and the Colour Gradient Model. 

\subsection{The Shan-Chen pseudo potential model}
  
We start with the standard LBM with Bhatnagar-Gross-Krook (BGK) collision term. In the lattice Boltzmann method, fictional particle groups on lattice nodes with discrete velocity are used to describe fluids. A distribution function $f_i(x)-f(x,e_i)$ is used to describe the occupation of each lattice site. Each lattice velocity $e_i$ on each lattice node $x$ has a distribution function $f_i(x)$. 




 Descretization of space, velocity and time are carried out in LBM. And this procedure greatly simplified the original Boltzmann equation. The location of particles in the space is restrained on the nodes of lattice grid, and the particle velocity is simplified into a very limited number of lattice velocities. We take a 2-D model as an example. This model is well known and widely used in the application. It is two dimensional and contain 9 velocities with the name D2Q9. This model is proposed by Qian et al \cite{qian1992}. In LBM, we assume all the particles has the same uniform mass (normally 1 mass unit is taken for simplicity). And the lattice unit($lu$) and time steps ($ts$) is important length and time unit in LBM. We only discuss with uniform mesh in this chapter ($\Delta x= \Delta y$).The macroscopic hydrodynamic quantities is calculated as:
	\begin{eqnarray}
	&\rho(\mathbf{x},t)=\sum_{i} f_i(\mathbf{x},t)& \\
	&\rho \mathbf{u}(\mathbf{x},t)=\sum_{i} \mathbf{c}_i f_i(\mathbf{x},t)&
	\end{eqnarray}


We use a D2Q9 lattice model \cite{qian1992} with nine velocities on a uniform squared two dimensional lattices in the simulation. This model is widely used in 2D lattice Boltzmann simulations. The lattice velocity could be written as:

\begin{equation}
\mathbf{e}=e\begin{bmatrix}
0 & 1 & 0 &-1 &0 &1& -1& -1& 1 \\
0 &0 &1 &0 &-1& 1& 1& -1& -1 
\end{bmatrix}
\end{equation}

Where $e=\Delta x / \Delta y$ is the local lattice speed, and have the relation with sound speed as: $c_s=\frac{c}{\sqrt{3}}$



The evolution of the distribution functions is described by the BGK collision terms:

\begin{equation}
f_i(\mathbf{x}+\mathbf{e}_i \Delta t,t+\Delta t)=f_i(\mathbf{x},t)-\frac{f_i(\mathbf{x},t)-f^{eq}_i(\mathbf{x},t)}{\tau}
\end{equation} 
The equilibrium distribution function is defined as:

\begin{equation}
f^{eq}_i(x)=w_i\rho (x)[1+3\frac{\mathbf{e}_i \cdot \mathbf{u}}{e^2}+\frac{9(\mathbf{e}_i \cdot \mathbf{u})^2}{2e^4}-\frac{3\mathbf{u}^2}{c^2}]
\end{equation}

Where the weight coefficients for D2Q9 model are:
\begin{equation}
w_i=\begin{cases}
4/9 & i=0 \\
1/9 & i=1..4 \\
1/36 & i=5..8
\end{cases}
\end{equation}

Macroscopic transport equations for mass, momentum and energy can be derived from Boltzmann equations using Chapman-Enskog expansion. \cite{ekexpansion} Kinematic viscosity $\nu$ in the D2Q9 model is obtained as:
\begin{equation}
\nu=c^2_s(\tau-\frac{1}{2})\Delta t
\end{equation}

An interaction force term between particles is used to describe the interparticle forces microscopically:

\begin{equation}
F(\mathbf{x},t)=-G\psi_{\sigma} (\mathbf{x},t)\sum_{i=1}^{8}w_i \psi_{\sigma} (\mathbf{x}+\mathbf{e}_i\Delta t,t)\mathbf{e}_i
\label{Shanchenforce}
\end{equation}

$\sigma$ and $\bar{\sigma}$ denotes the different fluid component. $G$ is a parameter determine the interaction strength $F(\mathbf{x},t)$ between neighboring particles. It also determine if the interaction is attraction or compulsion. To simulate a immiscible fluids system, $G$ should be kept a positive value so that a force will be generated to separate the fluids away from the interface. $\psi_{\sigma}$ is the interaction potential which is taken as component density : $\psi_{\sigma}=\rho_{\sigma}$. \\

\cite{shan2006} studied the spurious currents phenomenon of the Shan-Chen model. He conclude that the spurious currents phenomenons are due to the insufficient isotropy of the discrete gradient operator. A new method of obtaining the discrete gradient operator with higher order isotropy is given. However, more adjacent nodes are needed to calculate the gradient operator. It will increase the amount of massages that needed to be exchanged in the parallel computing. This action reduces the level of inherent parallelism of the lattice Boltzmann method.

\subsection{The Free Energy Model}
\cite{freeenergy} developed a thermodynamically consistent
binary fluid LB model by introducing an equilibrium state associated
with a free energy functional, corresponding pressure tensors and
chemical potentials. Correct choice of collision rules ensure that
the system evolutes toward to minimization of the free energy
functional.

The thermodynamic properties of a binary fluid system can be
described by a Laudau free energy functional.

\begin{equation}
\Psi=\int_V (\psi_b+\frac{\kappa}{2}(\partial_{\alpha} \phi)^2)dv
+\int_S \psi_s ds
\end{equation}

Where $\psi_b$ is the bulk free energy density and has a form:
\begin{equation}
\psi_b=\frac{c^2}{3} \rho ln\rho +A (-\frac{1}{2}\phi^2+
\frac{1}{4}\phi^4)
\end{equation}

$\phi$ is the order parameter represents the concentration of
components with the definition as:

\begin{equation}
\phi=\frac{\rho_a-\rho_b}{\rho_a+\rho_b}
\end{equation}

The hydrodynamics and thermodynamics of the binary fluids are
described by the Navier-Stokes equations and a convection-diffusion
equation:

\begin{equation}
\frac{\partial \phi}{\partial t}+\nabla (\phi u)= \mu \nabla ^2 \mu
\end{equation}

where the parameter $\mu$ determines the diffusivity of the binary fluids system.

A new distribution function $g_i(x,t)$ is introduced to describe the
concentration $\phi=\sum\limits_i g_i$ and is related to the
convection and diffusion. The distribution function $f_i(x,t)$ is
related to the fluid density and momentum as usual.

The time evolution equations using MRT scheme can be described as:
\begin{eqnarray}
\text{Collision Step:} & f'_i(x,t)=f_i(x,t)-M^{-1}SM(f_i-f^{eq}_i) & \\
& g'_i(x,t) =g _i^{eq}(x,t)& \\
\text{Streaming Step:} & f_i (x+e_i \Delta t,t+\Delta t) =f'_i(x,t)&
\\
&g_i (x+e_i \Delta t,t+\Delta t) =g'_i(x,t)&
\end{eqnarray}

Appropriate choice of equilibrium distribution functions can
reproduce the macroscopic equations in continuum limit. 

\subsection{The Color Gradient lattice Boltzmann Model}
An immiscible fluid model developed from Lattice Gas Cellular
Automata was introduced by \cite{color1991} The
particles in this model are colored either red or blue, therefore it
is normally called color gradient method. The surface tension is
introduced by adding a perturbation to the collision operator while
keep the adherence to the Navier-Stokes equations in homogeneous
regions. A recoloring step is involved after surface tension
perturbation calculation in order to achieve zero diffusivity of one
color into the other.

We use $f_i^r$, $f_i^b$ and $f_i$ denote the distribution functions
of a red fluid, a blue fluid and their combination respectively.A perturbation is computed to generate the a surface tension.The
surface tension can be expressed as a local anisotropy in the
pressure: the pressure measured normal to the surface is larger than
that tangential to the surface.The pressure in a LBM is proportional to the density, so surface
tension can be generated by placing preferentially particles in
directions normal to the interface rather than tangential. Mass and
momentum should be conserved.A color gradient G is defined as:

\begin{equation}
G(x,t)=\sum\limits_i e_i \sum\limits_j(f_j^r(x+e_i\Delta
t,t)-f_j^b(x+e_i\Delta t,t))
\end{equation}

Perturbation of the populations are:
\begin{equation}
f''_i(x,t) = f'_i(x,t) +A|G(x,t)|cos(2\theta)
\end{equation}


A redistribution of color enforce particles to move towards to the
regions occupied by particles with same color. This recoloring step
enable us to achieve the separation of the fluids. It is carried out
by the following maximization problem \cite{color1991}. Because the particles are forced to stay with the particles having the same colour, the Color Gradient model is only able to simulate complete immiscible fluids. 

\begin{equation}
W(f_i^{r''},f_i^{b''})=\max\limits_{f_i^{r''},f_i^{b''}}[\sum\limits_i
(f''_{r_i}-f''_{b_i})e_i]
\end{equation}

Subject to constrains:

\begin{eqnarray}
&\rho''_r=\sum\limits_i f_i^{r''}=\rho_r&\\
&f_i^{r''}+f_i^{b''}=f_i&
\end{eqnarray}

A general color gradient lattice Boltzmann model can be summarized
as:

\begin{itemize}
\item Single phase collision.
    %\begin{equation}
    %f_i=f_i^r+f_i^b
    %\end{equation}
\item Add a surface tension perturbation to $f'_i$ obtaining $f''_i$
\item Recoloring
\item Steaming
\end{itemize}

In both the Free Energy model and the Color Gradient model, only the value of distribution functions of nearest neighbouring nodes are needed, which keep the inherent parallelism of the lattice Boltzmann method.



\section{Comparison study}
We use three multi-component LB models for two problems : Poiseuille flow for binary fluids with viscosity contrast and capillary fingering. The numerical results are presented and compared with theoretical predictions. The thickness of the interface, limitation of the model and numerical stability are also discussed. A summary is given to show the capability and limitation of every model. 
 
\subsection{Poisseuille flow simulation for a immiscible binary fluid system with viscosity contrast}
Simulations of two fluid having a viscosity contrast in a channel is carried out and compared with theoretical predictions.(Figure\ref{CGFE_v}).The domain is periodic in y direction and bounded with non slip walls. The system is initialized with substance 0 in the middle and substance 1 on both sides near the walls. Initial density value 1 is applied to substance 0 and 1. The force acceleration applied to each substances is $10^{-6}$.\\

Figure(\ref{sc_f}) shows the velocity profile for the Shan-Chen model with viscosity ratio 12.  The simulation result does not match the analytical solution. Some noise is founded near the interface. To explore the reason of low agreement with the theoretical prediction, the density profile of two substances is shown as well. In the density profile, we can find big diffusion between two substances. The concentration of substance 0 decreased in the middle and increased in the domains that near the walls. Our previous study showed that good agreements with theoretical predictions can be obtained with fluids with viscosity ratio less than 10 \cite{chin_2002}. The flow of binary fluids with viscosity ratio larger than 10 cannot be calculated accurately, and the numerical instability emerge.  

\begin{figure}[h]
\begin{minipage}[c]{0.45\linewidth}
\centering
\includegraphics[width=\textwidth]{Shan-Chen_velocity.eps}
\end{minipage}
\hfill
\begin{minipage}[c]{0.45\linewidth}
\centering
\includegraphics[width=\textwidth]{Shan-Chen_density.eps}
\end{minipage}
\caption{Left: Simulated velocity of the Shan-Chen model; Right: Simulated density of the Shan-Chen model. Viscosity ratio between two substance is 10. The initial densities are both set as 1}\label{sc_f}
\end{figure}


\begin{figure}[h]
\begin{minipage}[c]{0.45\linewidth}
\centering
\includegraphics[width=\textwidth]{CG_poisseuille.eps}
%\caption{Simulated velocity of Color Gradient Model \label{CG_v}}
\end{minipage}
\hfill
\begin{minipage}[c]{0.45\linewidth}
\centering
\includegraphics[width=\textwidth]{FE_poisseuille.eps}
%\caption{Simulated velocity of Free Energy Model \label{FE_v}}
\end{minipage}
%\hfill
%\begin{minipage}[c]{0.28\linewidth}
%\centering
%\includegraphics[width=\textwidth]{FE_poisseuille.eps}
%\\hfillcaption{Simulated velocity of Free Energy Model \label{FE_v}}
%\end{minipage}
\caption{Simulated velocity of the Color Gradient Model (left) and the Free Energy Model (right) with viscosity ratio 100} \label{CGFE_v}
\end{figure}

The simulation results of both the Free Energy Model and the Color Gradient Model 
give excellent agreements with the analytical solution. The
concentration profile is also as expected.  Two immiscible fluid are
constrains in their respective areas and no unexpected diffusion is
discovered. We should notice that the interface thickness of the Free Energy Model is around $6$ lattice
units and the interface location is found at $x=38$ instead of their
initial position of $x=35$, a interface movement of 3 lattice is discovered from the result. The reason of this interface shifting is not clear and need further studies. The interface thickness of the Color Gradient Model is around 4, which is smaller than that of the Free Energy Model. A interface shifting is observed in color gradient model in
around 1.5 lattice units which is significantly smaller than the
free energy model.The recoloring algorithm separates two immiscible fluids very well
and gives nearly 0 diffusivity as we expected.\\

The thickness of the interface for different multi-component models varies. The thicknesses of the interface are studied in the Poiseuille flow simulations with a viscosity ratio of 10 and a surface tension of 0.05. The Shan-Chen model generate a thickness of 6 lattice unites, while the Free Energy model and the Color Gradient model gives 4 and 2-3 lattice units respectively.   

\subsection{Capillary fingering simulation}
Capillary fingering is a kind of hydrodynamic instability which exists in various displacements during oil/gas production.
 One fluids is pushed by another one, with different viscosities, along a
channel with non-slip wall.A growing finger of the driving fluid will be produced if the capillary
number $Ca$ is big enough.


\begin{equation}
Ca=\frac{u_t \rho \nu_2}{\sigma}
\end{equation}
where $u_t$ is the velocity of the tip of the finger, $\nu_2$ is the
viscosity of driving fluid. $\sigma$ is the surface tension.\\


Chin et al. studied the viscous fingering using the Shan-Chen model\ref{chin_2002}. Although fingering structures were observed in the simulation, the development of a finger is more clear in simulations at higher surface tensions, which is counterintuitive. We normally expect that increased surface tension is not favourable to the generation of fingerings. The author is not clear if this structure was produced due to viscous fingering or due to other effects.Numerical instability and high diffusivity was also founded in the simulation when the viscosity ratio is larger than 7. \\



We use the Color Gradient Model and the Free Energy Model to study the viscous fingering. A domain with grid sizes of 512x32 is used in the simulation. The first half of the domain contains substance 0 and the rest part is occupied by phase 1. Periodic boundary condition is used in X direction and bounded with non slip boundaries. A pressure gradient is imposed by applying a body force in X direction. Figure (\ref{fingering2}) shows the evolutions of fingers for binary fluids with a viscosity ratio of 20 and a tip velocity of 0.01. It is simulated by the Free Energy Model. From up to the bottom, they are finger evolutions for surfance tension 0.0678,0.0389,0.01985,0.00992 respectively. No finger will be produced if the surface tension is high. When the surface tension decrease, fingers are observed. However, the smaller surface tension, the less stable interfaces were produced. Figure (\ref{fingering3}) shows the evolution of the fingerings of the same binary fluids with a larger tip velocity of 0.05 that is calculated by the Color Gradient Model. Stable fingerings are observed in the results. \cite{halpern_1994} studied the fingering phenomenon in a channel.  The width of the fingers produced by different capillary number $Ca$ are measured.Their results are ploted in the figure (\ref{fingering1})with a solid line. A related width value with the width of the channel was used. Our results from the Free Energy Model and the Color Gradient Model are shown in the same figure with triangle points and star points respectively. 

\begin{figure}[h]
\begin{center}
\includegraphics[width=1.5in]{Capillary_Fingering1.eps}
\end{center}
\caption{From top to bottom: Fingering evolution for the Free Energy Model with surface tensions of 0.0678,0.0389,0.01985,0.00992 at a time interval of 1000 time steps. Viscosity ratio is 10; the tip velocity is 0.01}\label{fingering2}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[width=1.5in]{bw1e-4_2e-4_6e-4_1e-3.eps}
\end{center}
\caption{From top to bottom: Fingering evolution for the Color Gradient Model with surface tensions of 0.0389,0.01985,0.00992,0.00496 at a time interval of 1000 time steps. Viscosity ratio is 20; the tip velocity is 0.05}\label{fingering3}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[width=2.5in]{fingering.eps}
\end{center}
\caption{Fingering with as a function of Capillary number, the results from the Free Energy model is represented with triangle points, the results from the Color Gradient Model is represented with star points}\label{fingering1}
\end{figure}

Good agreements are achieved although some small discrepancies are found. These discrepancies might come from the boundary conditions. The Free Energy Model gives better agreements in low Capillary number simulations and the Color Gradient Models did better in high Capillary number simulations. These numerical examples shows that, fingering phenomena can be simulated properly by both the Free Energy Model and the Color Gradient Model captured. It is worth mentioning that, the numerical stability for both the Free Energy Model and the Color Gradient Model is similar, the maximum viscosity ratio for dynamic interfaces simulations is around 20. \\

\section{Summary}
The feathers and limitations of the Shan-Chen model, the Color Gradient Model and The Free Energy Model  are investigated on seven aspects: The maximum density ratio of binary fluids that can be simulated, the maximum viscosity ratio of binary fluids that can be simulated, interface thickness in the poiseuille simulation, interface shifting distance in the poiseuille simulation, level of inherent parallelism, ability of simulating capillary fingering and diffusion property. \\

The Shan-Chen model can simulate high density ratio which is up to 800 binary fluids with same viscosity. Viscosity contrast will lead to high diffusive interface in the Shan-Chen model and therefore affect significantly the numerical stability. The Free Energy model and the Color Gradient model has similar capability on this point: they can simulate binary fluids having same density with viscosity contrast. In our study, we found that, the maximum viscosity ratio of these two models depends on the type of problem. If the interface is static, the fluids with viscosity ratio up to 120 can be simulated, however, for dynamic interface simulation, this value decrease to 20. 
\begin{center}
%\small%\addtolength{\tabcolsep}{-15pt}
\begin{tabular}{|>{\columncolor[rgb]{0.8,0.8,0.8}}c|ccc|}

\hline
% X/Y &  \multicolumn{4}{>{\columncolor[rgb]{0.8,0.8,0.8}}c|}{Poblacin } \\ 
 \rowcolor[rgb]{0.8,0.8,0.8}  & \multicolumn{1}{p{3cm}}{\centering {\scriptsize The Shan-Chen Pseudo Potential Model}} & \multicolumn{1}{p{3cm}}{\centering \scriptsize The Free Energy Model} &\multicolumn{1}{p{3cm}|}{\centering  \scriptsize The Color Gradient Model} \\ 
\hline 
 \multicolumn{1}{>{\columncolor[rgb]{0.8,0.8,0.8}}p{3cm}|}{\centering \scriptsize Maximum density ratio} & 800 &  1 &  1  \\ \hline

 \multicolumn{1}{>{\columncolor[rgb]{0.8,0.8,0.8}}p{3cm}|}{\centering \scriptsize Maximum viscosity ratio} &  5 & \multicolumn{1}{p{3cm}}{\centering \scriptsize 120 (static interfaces) \\ 20 (dynamic interfaces) } & \multicolumn{1}{p{3cm}|}{\centering \scriptsize 120 (static interfaces) \\ 20 (dynamic interfaces) }  \\ \hline 
 
 \multicolumn{1}{>{\columncolor[rgb]{0.8,0.8,0.8}}p{3cm}|}{\centering \scriptsize Interface thickness (Poiseuille flow simulation)}  & 10 $lu$ & 5 $lu$ & 2-3 $lu$
\\ \hline

 \multicolumn{1}{>{\columncolor[rgb]{0.8,0.8,0.8}}p{3cm}|}{\centering \scriptsize Interface Shifting distance (Poiseuille flow)}  & 10 $lu$ & 5 $lu$ & 2-3 $lu$
\\ \hline

\multicolumn{1}{>{\columncolor[rgb]{0.8,0.8,0.8}}p{3cm}|}{\centering \scriptsize Ability of simulating capillary fingering}  & No & Yes & Yes
\\ \hline

\multicolumn{1}{>{\columncolor[rgb]{0.8,0.8,0.8}}p{3cm}|}{\centering \scriptsize Inherent Parallelism}  & Low  & High & High
\\ \hline

\multicolumn{1}{>{\columncolor[rgb]{0.8,0.8,0.8}}p{3cm}|}{\centering \scriptsize Diffusion}  & {\scriptsize High diffusivity interface}  & \multicolumn{1}{p{3cm}}{\centering \scriptsize Diffusivity can be controlled by a single parameter}  & \multicolumn{1}{p{3cm}|}{\centering \scriptsize Complete immiscible, no diffusion on the surface} 
\\ \hline

\end{tabular}\label{comp}
\end{center}


%% The Appendices part is started with the command \appendix;
%% appendix sections are then done as normal sections
%% \appendix

%% \section{}
%% \label{}

%% References
%%
%% Following citation commands can be used in the body text:
%%
%%  \citet{key}  ==>>  Jones et al. (1990)
%%  \citep{key}  ==>>  (Jones et al., 1990)
%%
%% Multiple citations as normal:
%% \citep{key1,key2}         ==>> (Jones et al., 1990; Smith, 1989)
%%                            or  (Jones et al., 1990, 1991)
%%                            or  (Jones et al., 1990a,b)
%% \cite{key} is the equivalent of \citet{key} in author-year mode
%%
%% Full author lists may be forced with \citet* or \citep*, e.g.
%%   \citep*{key}            ==>> (Jones, Baker, and Williams, 1990)
%%
%% Optional notes as:
%%   \citep[chap. 2]{key}    ==>> (Jones et al., 1990, chap. 2)
%%   \citep[e.g.,][]{key}    ==>> (e.g., Jones et al., 1990)
%%   \citep[see][pg. 34]{key}==>> (see Jones et al., 1990, pg. 34)
%%  (Note: in standard LaTeX, only one note is allowed, after the ref.
%%   Here, one note is like the standard, two make pre- and post-notes.)
%%
%%   \citealt{key}          ==>> Jones et al. 1990
%%   \citealt*{key}         ==>> Jones, Baker, and Williams 1990
%%   \citealp{key}          ==>> Jones et al., 1990
%%   \citealp*{key}         ==>> Jones, Baker, and Williams, 1990
%%
%% Additional citation possibilities
%%   \citeauthor{key}       ==>> Jones et al.
%%   \citeauthor*{key}      ==>> Jones, Baker, and Williams
%%   \citeyear{key}         ==>> 1990
%%   \citeyearpar{key}      ==>> (1990)
%%   \citetext{priv. comm.} ==>> (priv. comm.)
%%   \citenum{key}          ==>> 11 [non-superscripted]
%% Note: full author lists depends on whether the bib style supports them;
%%       if not, the abbreviated list is printed even when full requested.
%%
%% For names like della Robbia at the start of a sentence, use
%%   \Citet{dRob98}         ==>> Della Robbia (1998)
%%   \Citep{dRob98}         ==>> (Della Robbia, 1998)
%%   \Citeauthor{dRob98}    ==>> Della Robbia


%% References with bibTeX database:

\bibliographystyle{elsarticle-harv}
%\bibliographystyle{elsarticle-num}
\bibliography{ref}

%% Authors are advised to submit their bibtex database files. They are
%% requested to list a bibtex style file in the manuscript if they do
%% not want to use elsarticle-harv.bst.

%% References without bibTeX database:

% \begin{thebibliography}{00}

%% \bibitem must have one of the following forms:
%%   \bibitem[Jones et al.(1990)]{key}...
%%   \bibitem[Jones et al.(1990)Jones, Baker, and Williams]{key}...
%%   \bibitem[Jones et al., 1990]{key}...
%%   \bibitem[\protect\citeauthoryear{Jones, Baker, and Williams}{Jones
%%       et al.}{1990}]{key}...
%%   \bibitem[\protect\citeauthoryear{Jones et al.}{1990}]{key}...
%%   \bibitem[\protect\astroncite{Jones et al.}{1990}]{key}...
%%   \bibitem[\protect\citename{Jones et al., }1990]{key}...
%%   \harvarditem[Jones et al.]{Jones, Baker, and Williams}{1990}{key}...
%%

% \bibitem[ ()]{}

% \end{thebibliography}

\end{document}

%%
%% End of file `elsarticle-template-harv.tex'.
